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Abstract 

We recommended consequent discrete combinatorial research in math- 
ematical physics 6 . Here we show an example how discretization 
of partial differential equations can be done and that quickly unex- 
pected new findings can result from research in this up to now un- 
explored area. We transformed the vacuum Maxwell equations into 
finite-difference equations, provided simple initial conditions and stud- 
ied the development of the electromagnetic fields using special software 
[7j jSJ. The development is wave-like as expected. But it is not trivial, 
the wave maxima have different heights. If all (by definition minimal) 
finite differences of the location coordinates are multiplied by num- 
bers (coupling factors) whose squares are equal to the fine structure 
constant, we noticed: 

1. The first two wave maxima have nearly the same height. Of course 
this can be also coincidental. 

2. The following maxima are at first slightly decreasing and then, 
beginning with the 6th maximum, exponentially increasing. 



PACS: 42.25.Bs, 02.70.Bf, 02. 10. Ox 
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1 Introduction 



We have shown jH] that only an a priori finite mathematical model can have 
an exact equivalent in physical reality. This means that it implies an only 
finite number of basic arithmetic operations on an a priori finite numerical 
space 1 which can be represented without using irrational numbers. 

Up to now there is not much experience in this area: Important physical 
equations are defined on continuous (a priori infinite) sets and often written 
as partial differential equations. If we want to find the natural finite ba- 
sis of them, first we have to replace differential calculus by finite-difference 
calculus. This can soon lead to difficult combinatorics, especially in case of 
interactions across several dimensions. But increasing performance of com- 
puters offers new possibilities. The mentioned numerical space can be rep- 
resented by finite dimensional numerical lattices (finite dimensional point 
lattices with numbers resp. quantities assigned to every point) which can 
be handled adequately by a computer. So we developed as help software for 
handling numerical lattices and for studying the results of numerical algo- 
rithms on them. All lattice points are addressed by integer coordinates and 
the quantities assigned to the points can be complex 2 . We used the software 
to study the vacuum Maxwell equations. 

We converted these into finite-difference equations as directly as possible be- 
cause we wanted to operate on an elementary level (a priori finite). So we 
had to avoid analytic methods (like usage of infinite power series, trigonomet- 
ric functions, complex and matrix exponential functions, Fourier transforms) 
which are often used to force an easy to survey, stable development. 

2 Material and Method 
2.1 Software 

The in C++ written open-ended software [7] [S] is designed for editing of 
numerical lattices and for generation of statistics like sums over the absolute 

1 Nevertheless both can increase without boundary when time increases without bound- 
ary (infinite potential). 

2 Both complex rational and complex floating-point numbers are supported. Subse- 
quently by the term lattice always this kind of numerical lattice is meant. 



3 



or squared quantities of subspaces, with graphical representation. Using the 
ASCII file format of saved data as interface or writing own code the user 
can check the results of own algorithms. Additionally it is possible to define 
algorithms without writing code. This is done by usage of a table with sev- 
eral entries. Every entry defines a coupling between two lattices in form of 
an addition from the first (the source lattice) to the second (the destination 
lattice). The strength of the coupling is determined by the coupling factor p: 
During every iteration of the algorithm all non-zero quantities of the source 
lattice with maximal iteration number 3 are added to shifted locations within 
the destination lattice after multiplication by p. 

Source and destination lattices, the multidimensional relative shift resp. off- 
set are selectable using integer numbers, the coupling factor p can be a com- 
plex number, if desired. The user determines how many iterations of all 
additions (defined by the table) are done. 

The motivation for this scheme arises from the fact that it is simple and nevertheless very flexible. Among 
others it can be used to study random walks like those discussed in UJ. Note that p also can be a complex 
probability amplitude and is not restricted to the interval [0, 1]. 

If the number of iterations and additions per iteration is great, the lattice(s) 
may occupy too much memory of the computer. Therefore it is possible 
to activate automatic deletion of the lattice points with smallest absolute 
quantities if these would occupy too much memory. Here we used this after 
we have checked that up to 150 iterations in our application it has not affected 
the 12 most significant digits of the evaluated numerical data. 

2.2 The Maxwell equations 

Using MKS units, the Maxwell equations in linear and isotropic media are 

f. CD 
0, (2) 

V x B , (3) 

-V x E , (4) 

where E is the electric field, B is the magnetic field, p is the charge density, J 
is the vector current density, e is the permittivity and p is the permeability. 

3 The results of every iteration are indexed by an increasing iteration number n. This 
can be used to analyze the temporal development. 



V-E 



V B 

dE 

dB 

~dt 
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In the vacuum which is free of current and electric charge (J = and p = 0) 
they reduce to the vacuum Maxwell equations in which besides (J2J) also ((H) 
vanishes and (j3J), (jlj) become 

— VxB = c 2 VxB, 

/^oeo 

—V x E , 

where eo is the permittivity of free space, /io is the permeability of free space 
and c is the speed of light in the vacuum. Using the definition 

B := cB 

we obtain 

c V x B , (5) 

-c V x E . (6) 

The proportionality factor c still has the dimension ms" 1 . We will see (cf. 
12. 4j) that during discretization automatically the question about the natural 
dimensionless relationship arises. 



<9E 
~dt 
dB 
~dt 



dE 
~dt 
dB 
~dt 



2.3 Discretization 

Let x,y,z denote Cartesian coordinates and let E x ,E y ,E z denote the com- 
ponents of E = E(t, x, y, z) and B x , B y , B z those of B = B(t, x, y, z). With 
this © and (0) are equivalent to 

, (7) 



dE x _ J dB z 8B y \ dE y _ J dB x dB z \ dE z _ J dB y dB x 
dt \ dy dz J ' dt \ dz dx I ' dt \ dx dy 



dB x fdE z dE y \ dB y ( 8E X dE z \ dB z ( dE v dE., 



dt V dy dz J dt V dz dx J dt V dx dy 



(8) 



A careful consideration, however, must take into account that the differences 
in time and location cannot be arbitrarily small. 



5 



2.3.1 (Half) minimal differences as units 

We know that there are non- vanishing minimal absolute differences dt in time 
and ds in location 4 so that absolute differences of t below dt and of x, y, z 
below ds are not existing in physical reality. Therefore it is appropriate to 
base the scaling on dt and ds. We will use dt as unit of time and (because 
we want to form minimal symmetric differences in location) ds/2 as unit of 
location. 



2.3.2 Integer representation 

Define 



* . 2a; A 2y A 2z 

1 := 17 > x := T> V := T> z := T and 
dt ds ds ds 



(9) 



E(t, x, y, z) := E(t, x, y, z), B(t, x, 5) := B(t, x, y, z) . (10) 

The coordinates t, x, y, z are all dimensionless having only integer differences. 
Now it is possible to form minimal symmetric finite differences using only 

integer coordinates. For v G {x, y, z} let E v , B v , E v and B v denote the 

components of E, B, E and B respectively. Define 

E v (i,x,y, z) = E v (t + dt,x, y, z) - E v (t, x, y, z) , 
E v (i, x — l,y, z) — E v (t, x + ds/2, y, z) — E v (t, x — ds/2, y, z) , 
E v (i,x,y— 1, z) = E v (t,x,y + ds/2, z) - E v (t,x,y - ds/2, z) , 
E v (i,x,y, z - 1) = E v (t,x,y, z + ds/2) - E v (t,x,y, z - ds/2) , 

B v (t, x, y, z) = B v (t + dt,x, y, z) - B v (t, x, y, z) , 
B v (t, x — l,y,z) = B v (t, x + ds/2, y, z) - B v (t,x- ds/2,y, z) , 
B v (t, x, y — 1, z) = B v (t, x, y + ds/2, z) — B v (t, x, y — ds/2, z) , 
B v (t, x,y, z — 1) = B v (t, x,y,z + ds/2) - B v (t,x, y,z- ds/2) , 

and transformation of (J7|) and (jHJ) into finite-difference equations yields 



E v t 


= E v (i+l,x,y,z 


Eyx 


= E v (i,x + l,y,z 


E 


= E v (i,x,y + l,z 


E vz 


= E v (i,x,y,z+\ 


B v t 


= B v (i+l,x,y,z 


Bvx 


= B v (t,x+l,y,z 


Byy 


= B v (i,x,y + l,z 


B vz 


= B v (i, x,y, z + 1 



E x t 
~~dT 



B 



ds 



B 



ds 



Eyt_ 

dt 



I B xz 


B zx \ 


E zt | 


( B 


B X y \ 


I ds 


ds J 


' -dt = C \ 


, ds 


ds ) 



B x t _ I E zy E yz \ B yt _ / E xz E zx \ B ' zt _ I E yx E ' xy \ 
dt \ ds ds } ' dt \ ds ds I ' dt \ ds ds I 



4 They depend on further variables like energy per photon which influence both dt and 
ds (but not necessarily in all those cases the quotient dt/ds). 
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This would be equivalent to ((7J) and (JSJ) in the borderline case dt — > and 
ds — > 0. But due to 12. 3. li the equations (fTTj) and (|T2|) are closer to reality. 
Using 

dt 

p := c 



B X y 



ds 

as absolute coupling factor they become more compact: 

-Ezt = P [B zy - B yz ^j , E yt = p (j3 xz ~ B zx ^j , E ' zt = p (i? 
B x t = — j5 (-^22/ — A/z) 7 Byt — —p [E X z — E zx ^j , B zt = —p {Eyx — Exy 

2.4 The algorithm 

Written out explicitly from (|14|) and (|15j) follows 



E x (i+l,x,y,z) = E x (i,x,y,z) + 



E y (i+l,x,y,z) = E y (i,x,y,z) + 



E z (i + l,x,y,z) = E z (i,x,y,z) + 



B x (i+ l,x,y,z) = B x (i,x,y,z) 



B y (i+ l,x,y,z) = B y (i,x,y, 5) 



B z (i+ l,x,y,z) = B z (i,x,y,z) - 



p[B z (t,x,y+ 1,5 
p (B y (i,x,y, z + 1 



p[B x (t,x,y,z + 1 
p(B z (i,x + l,y,z 
p(B y (t,x + l,y,z 



P 



(b x {1 



x,y + l,z 



p(E z (i,x,y + 1,5 
p (^E y (i,x, y, 5 + 1 
p (^E x {i,x,y, z + 1 
p(^E z (i,x + l,y,z 
p(E y {i,x + \,y,z 
p(E x (i,x,y + l,z 



B z (i,x,y - 1,5 
B v (t,x,y,z- 1 
B x (i,x,y,z - 1 
B z (i,x- 1,2/, 5 
B y (i,x~ l,y,z 
B x (i,x,y- 1,5 

E z (i,x,y- 1,5 
E y (i,x,y,z - 1 
E x (i,x,y, 5-1 
E z (i,x- 1,2/, 2 
E y (i,x- l,y,z 
E x (i,x,y- 1,5 



(13) 

(14) 
(15) 



(16) 



For stepwise calculation of the temporal development we can, as shown in [7] , 
directly implement these equations in our software using 6 four- dimensional 
numerical lattices, 3 representing E, 3 representing B, and ±p as coupling 
factor p (cf. 12.1)1 . At this p is a free variable, so we have to search for 
reasonable values of it (i.e. automatically we have to search for the natural 
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dimensionless relationship, as mentioned in !2.2j) . The time steps are uniform 
because both laws (J5J) and (jUJ) are considered to be valid at every time, i.e. 
E and B are calculated simultaneously in (fT^j) . not alternating. Therefore 
this algorithm differs from leapfrog FDTD schemes [TU] . 

The above equations show that a linear superposition of the initial con- 
ditions leads to a linear superposition of the results. Therefore we chose 
for the purpose of clarity as simple as possible nontrivial (non-zero) initial 
conditions. As initial time we defined t = and set all initial quantities 
(of E x , E y , E z , B x , By, B z ) to zero on all locations except E x in the origin 
(0,0,0,0) which we set to one: E x (0, 0,0,0) := 1. Then we analyzed the 
temporal development of E x resulting from iterative calculation of and 
its dependence of p which by definition (fTHj) is a positive dimensionless num- 
ber. 

3 Results 

As expected a wavelike development results. But it is not trivial, the heights 
of the wave maxima are not constant. In case of p > 1/8 they are increasing 
from the beginning, which is not realistic. Therefore we focused our interest 
on smaller values of p. 

3.1 Preliminary note 

Careful evaluation is necessary because (j!6|) bases on a simplified model, 
particularly it is derived from the idealized vacuum situation. Even small 
errors in can become relevant, if many iterations are done. Therefore 
the calculation of the first wave is more reliable than calculations of other 
waves and a detailed graphical representation of the first wave is reasonable. 

3.2 Short term development: The first wave 

We are searching for realistic values of p. If these should cause an as stable 
as possible wave pattern, the wave maxima should have nearly equal height. 
Especially the first two maxima (the maxima of the first wave which we can 
calculate with the best reliability) should have similar heights. Fig. ^ |2] and 
EH show the first wave for three different values of p. 
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Figure 1: The first wave in case of p=l/16. Abscissa: t, 
Ordinate: E x (t, 0, 0, 0). 




Figure 2: The first wave in case of p — \fu = 0.085424542921. Abscissa: t, 
Ordinate: E x (t,0,0,0). 




Figure 3: The first wave in case of p—1/8. Abscissa: t, 
Ordinate: £^(4,0,0,0). 
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We see in Fig. [T] (resp. EJ) that in case of p = 1/16 (resp. p = 1/8) the 
second maximum is considerably smaller (resp. greater) than the first. More 
interesting is Fig. El There 

p = ^ = 0.085424542921 = 1/11.7062376432 

in which a = 1/137.03599976 is the fine structure constant 5 as listened in 

The calculated quantities near the maximum at t = 28 are 

4(26,0,0,0) = 0.877636902081288, 

4.(27,0,0,0) = 0.950716819197347, 

4(28,0,0,0) = 0.98752930561647, 

4(29,0,0,0) = 0.986271782354442, 

4(30,0,0,0) = 0.947405056005354. 

So in case of p = \/a the second maximum deviates less than 2% from the 
first maximum which lies in the start and has the initial quantity 1. 6 

3.3 Long term development 

We see in Fig. 0J |3] and El that E x (t, 0, 0, 0) begins to diverge considerably 
from the 6th maximum on. Fig. [7|and|B]show the development of E x (t, 0, 0, 0) 
together with B z (t, 0,1,0) and Fig. El the local distribution of E x after 150 
iterations of ()16|). 



5 It can be derived very precisely from basic experiments, e.g. from measurements of 
the Quantum Hall Effect [TT], 

6 In case of p = 0.08672495 the maximum also is at t = 28 and differs less than 1 ppm 
from 1. But due to there is no good foundation for such fine interpolation. 
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Figure 4: Long term development in case of p — 1/16. Abscissa: t, 
Ordinate: E x (t,0,0,0). 
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Figure 5: Long term development in case of p — y/a. Abscissa: t, 
Ordinate: E x (t, 0,0,0). 
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Figure 6: Long term development in case of p—l/S. Abscissa: t, 
Ordinate: E x (t,0,0,0). 
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Figure 7: E x together with B z in case of p 
Ordinate: E x (t, 0, 0, 0) and B z (t, 0, 1, 0). 



y/a. Abscissa: t, 




Figure 8: Relation between the fields shown in Fig. [7| {p = y/a, t 6 [0, 160]). 
Abscissa: B z (t, 0, 1, 0), Ordinate: E x (t, 0, 0, 0). 




Figure 9: Local distribution of E x after 150 iterations of (fT7)|) in case of 
p— y/a. Abscissa: x, Ordinate: E x (150,x,0,0). 
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4 Comment 



It was not surprising that after many iterations of (fTfijl a wave packet results 
which is concentrated round the origin (Fig. 05)). As long as there is one 
and the same p for all directions, no direction is preferred. However, it was 
surprising for us that beginning with the 6th maximum there must be a 
dominating positive feedback from the surrounding which causes divergence 
of E x (t, 0, 0, 0) (and of all other field components in the origin). We noticed 
this also for p < 1/16. Further numerical tests indicated in case of p = ^Ja 
and t G [200, 1000] an exponential increase of the wave maxima with an 
average growth factor of about 1.04 per step t — > t + 1. Of course it is 
possible to force a stable development and conservation of energy by simple 
normalization, e.g. by division of the right sides of ()16|) by an appropriate 
function U{t) > 1, but we would like to have a good justification for doing 
that. Remembering 13.11 we guess that more advanced changes of (I16J) are 
necessary. Although the equations (fTtjj) are not exact, they can be suitable for 
approximative calculation of the initial discrete short term development (see 
13. 2|) because they are directly derived from the Maxwell equations. We have 
done this without any usage of our knowledge of the fine structure constant, 
and it is noteworthy that the first maxima of E x (t, 0, 0, 0) are nearly equal as 
shown in Fig. |21 if the in (J13|) defined absolute coupling factor p is the root 
of the fine structure constant. 



5 Conclusions 

The (vacuum) Maxwell equations cannot describe physical reality exactly, 
like all differential equations [5]. Straightforward conversion of the vacuum 
Maxwell equations into finite-difference equations leads to a system of equa- 
tions (HHJ) which is not exact, too. But it is unambiguous, provides additional 
combinatorial details and can be suitable for approximative calculation of the 
initial discrete short term development of the electromagnetic fields. The ab- 
solute coupling factor p is a new variable which automatically arises during 
formation of the finite differences. In case of p = ^J~a a short term wave-like 
development with nearly equal initial maxima results. If this is not coin- 
cidental, the definition of p (cf. ffT^j) together with 12.3. 1|) suggests a new 
interpretation of the fine structure constant a and ()16j) can serve as starting 
point for improvements and further combinatorial studies. 
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